%% This program solves for A,B, p1star and p2star, given certain inputs

function alphapstar = alphapstar1(alpha, gamma, meanp, beta, nonoilgdp, oilprod)

global pmuh pmul muh mul xmax expon Vaut

simequ0 = xmax';
options = optimset('TolFun',1e-6,'TolX',1e-6,'MaxFunEvals',10000,'MaxIter',25000,'Display','off');

% x(1) = p1star, x(2) = p2star, x(3) = capA, x(4) = capB
             
simequ  = @(x)  [(1/expon)*(nonoilgdp + oilprod*(alpha + gamma*(x(1) - meanp))).^expon + beta*x(3) + beta*pmuh*x(4) +  beta*pmuh*quad(@Ve_h,x(2),10000) + beta*pmul*quad(@Ve_l,x(1),10000) - ...  % V_h(p*) = V_e(p*,mu_l)
                                                                                                                            (1/expon)*(nonoilgdp + oilprod*x(1)).^expon + mul - beta*Vaut;
                 (1/expon)*(nonoilgdp + oilprod*(alpha + gamma*(x(2) - meanp))).^expon + beta*x(3) + beta*pmuh*x(4) +  beta*pmuh*quad(@Ve_h,x(2),10000) + beta*pmul*quad(@Ve_l,x(1),10000) - ...  % V_h(p*) = V_e(p*,mu_l)
                                                                                                                            (1/expon)*(nonoilgdp + oilprod*x(2)).^expon + muh - beta*Vaut;
                 x(3) - quad(@(z)VhInt(z,x(1),x(2),x(3),x(4),alpha,gamma),0,x(1));        % A = int(V_h)
                 x(4) - quad(@(z)VhInt(z,x(1),x(2),x(3),x(4),alpha,gamma),x(1),x(2))];    % B = int(V_h) 
             
xmax = fsolve(simequ,simequ0,options);                       

alphapstar = [xmax(1) xmax(2) xmax(3) xmax(4)];